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Abstract 


Finite  difference  schemes  do  not  perform  well  in  simulations  of  high  Reynolds  number  flow 
because  a  restrictive  number  of  grid  points  must  be  used  to  resolve  boundary  layers.  In  this 
report  we  present  a  grid-free  numerical  method  which  may  be  used  to  calculate  flows  around 
an  immersed  elastic  structure.  Numerical  results  are  presented  in  two  simple  cases:  flow 
past  a  circular  cylinder  and  flow  past  a  flat  plate.  Here  the  boundary  motion  is  specified  and 
is  not  time  dependent.  However,  this  algorithm  can  be  used  to  simulate  more  complicated, 
time  dependent  problems  in  biological  fluid  mechanics. 
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1.  Introduction 

This  project  concerns  the  computational  modeling  of  problems  in  biological  fluid  mechan¬ 
ics.  Typically,  these  involve  the  study  of  flow  about  highly  flexible,  moving  boundaries;  for 
instance,  distensible  vessels,  moving  heart  walls  and  valves.  The  fluid  along  with  the  bound¬ 
ary  constitues  a  coupled  mechanical  system  in  the  sense  that  the  motion  of  the  boundary  is 
determined  by  that  of  the  fluid,  but  at  the  same  time  the  boundary  exerts  force  on  the  fluid 
and  alters  its  motion.  With  the  advent  of  supercomputers,  many  of  these  problems  which 
were  previously  intractable,  can  now  be  successfully  analyzed. 

An  innovative  computational  approach,  the  immersed  boundary  technique ,  was  introduced 
by  Peskin  [16]  to  model  two-dimensional  blood  flow  in  the  left  heart.  Recently,  this  technique 
has  been  advanced  to  study  other  biofluiddynamic  problems  such  as  platelet  aggregation  [10] 
aquatic  animal  locomotion  [7] [9],  peristaltic  pumping  of  solid  particles  [8],  fluid  flow  in  the 
inner  ear  [3],  and  three-dimensional  blood  flow  in  the  heart  [17][18].  This  numerical  method 
solves  the  full  incompressible  Navier-Stokes  equations  in  a  domain  of  fluid  within  which 
a  massless,  neutrally-buoyant  elastic  boundary  undergoing  time  dependent  movements  is 
immersed. 

The  Navier-Stokes  equations  have  been  solved  using  Choinrs  iu:,U  difference  scheme  [5]. 
Theoretically,  this  scheme  imposes  no  restriction  on  the  Reynolds  number  of  the  flow  to  be 
modelled.  (The  Reynolds  number  is  a  nondimensional  quantity  which  measures  the  ratio  of 
inertial  forces  to  viscous  forces.)  However,  the  higher  the  Reynolds  number,  the  smaller  the 
grid  spacing  and  time  step  must  be  to  ensure  stability  of  the  calculations.  To  model  high 
Reynolds  number  flow  (as  in  the  heart),  the  amounts  of  computer  time  and  memory  needed 
become  prohibitive  and  impractical. 

Within  the  finite  difference  framework,  we  have  adapted  the  second  order  projection 
method  recently  introduced  by  Bell,  Colella,  and  Glaz  [2]  for  the  solution  to  the  Navier- 
Stokes  equations.  Our  results  indicate  that  flows  with  higher  Reynolds  numbers  can  be 
modelled  with  this  code  than  with  C'horin’s  scheme  using  the  same  number  of  grid  points 
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and  the  same  time  step.  We  have  tested  this  scheme  on  flow  past  a  tethered  cylinder  and 
have  exhibited  flow  separation  and  vortex  shedding.  The  cylinder  is  modelled  as  an  immersed 
boundary.  Samn  [20]  has  tested  the  projection  method  on  problems  with  a  non-zero  force 
density,  as  is  the  case  of  immersed  boundary  calculations,  in  a  class  of  problems  where  exact 
solutions  are  known.  For  moderately  high  Reynolds  numbers,  this  projection  method  is  quite 
promising. 

For  much  higher  Reynolds  numbers,  a  grid-free  vortex  method  combined  with  the  im¬ 
mersed  boundary  technique  may  be  more  appropriate.  An  effort  in  this  direction  was  made 
by  McCracken  and  Peskin  [15]  in  1980.  This  was  a  hybrid  method  which  represented  vor- 
ticity  both  on  a  grid  and  as  free-moving  vortex  elements.  In  light  of  the  advances  made  in 
vortex  methods  and  the  incredible  growth  of  computing  power  since  then,  a  totally  grid-free 
vortex  method  for  two-dimensional  flow  is  now  being  implemented  and  tested.  This  report 
will  describe  the  current  status  of  this  algorithm. 

In  section  2,  we  will  briefly  discuss  the  immersed  boundary  technique,  and  in  section  3  we 
will  describe  vortex  methods.  In  section  4  we  will  present  our  grid-free  vortex  method  and 
discuss  how  it  treats  an  immersed  boundary.  Examples  of  flow  past  a  cylinder  and  vortex 
roll-up  due  to  the  motion  of  a  flat  plate  will  be  presented  in  Section  5.  Finally  in  Section 
6  we  will  discuss  the  issues  which  need  to  be  addressed  to  make  this  method  robust  and 
competitive. 
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2..  Immersed  Boundary  Technique 

The  flow  of  a  viscous,  incompressible  fluid  is  governed  by  the  incompressible  Navier  Stokes 
equations: 


Vp  +  n  Au  +  F(x,  t)  (1) 

0 

Here  p  =  density,  p.  =  viscosity,  u  =  (u,  v )  =  velocity,  p  =  pressure,  and  F(x,  t)  is  the 
external  force  per  unit  volume  applied  to  the  fluid.  This  external  force  field  represents  the 
force  of  the  elastic  boundary  on  the  fluid.  It  is  a  delta-function  layer  of  force  which  is  zero 
away  from  the  immersed  boundary  X(s,t).  This  representation  will  be  the  basis  of  the 
computational  model. 

Since  the  immersed  boundary  is  taken  to  be  elastic  and  massless,  the  density  of  the 
boundary  force  f(s,f)  at  a  material  point  is  determined  by  the  boundary  configuration  at 
time  t.  We  will  discuss  the  derivation  of  the  boundary  force  below.  This  elastic  force 
is  transmitted  directly  to  the  fluid  because  the  boundary  is  massless.  This  follows  from 
Newton’s  second  law  for  force  balance  at  a  boundary  point:  the  force  of  the  fluid  on  the 
boundary  point  and  the  elastic  force  on  the  boundary  point  must  add  up  to  zero.  Therefore, 
the  force  of  the  boundary  point  on  the  fluid  is  equal  to  the  elastic  force  on  that  point.  The 
external  force  field  is  then: 


8u 

dt 


+  u  •  Vu 


V  •  u  = 


F(x,f)  =  J  f(s,t)8(x-X(s,t))ds  (2) 

Here  the  integration  is  over  the  immersed  boundary,  and  6  is  the  two-dimensional  delta- 
function. 

Since  the  fluid  is  viscous,  the  velocity  field  is  continuous  across  the  boundary.  This 
implies  that  the  velocity  of  a  material  point  of  the  immersed  boundary  must  be  equal  to  the 
velocity  of  the  fluid  at  that  point: 
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-  u(X(s,t),t)  =  Ju(x,t)S(x-X(s,t))dx  (3) 

Here  the  integration  is  over  the  entire  fluid  domain.  This  integral  representation  is  not 
singular  since  the  two-dimensional  delta-function  is  integrated  over  both  space  dimensions. 

The  crucial  feature  of  this  model  is  that  the  immersed  boundary  is  not  our  computational 
boundary  in  the  Navier-Stokes  solver  (whether  we  are  using  finite  differences  or  vortex  meth¬ 
ods).  It  is  a  singular  force  field  which  alters  the  driving  force  in  the  fluid  dynamic  equations. 
Algorithm 

Given  the  fluid  velocity  field  un  and  the  configuration  of  the  elastic  boundaries  X”  at  time 
step  n: 

1)  Calculate  the  force  density  fn  from  the  configuration  Xn. 

2)  Use  fn  to  determine  the  external  force  on  the  fluid:  Fn. 

3)  Solve  the  Navier-Stokes  equations  for  un+1. 

4)  Move  the  elastic  boundary  at  the  local  fluid  velocity  to  get  Xn+1. 

Within  the  finite  difference  framework,  step  (3)  has  been  solved  using  Chorin’s  finite 
difference  method.  Steps  (2)  and  (4)  involve  the  use  of  discrete  delta  functions  which  com¬ 
municate  information  between  the  grid  and  the  immersed  boundary  points  (which  do  not 
coincide  with  grid  points).  We  refer  the  reader  to  [16]  for  details. 

Derivation  of  force  density 

We  will  first  model  a  stationary  immersed  boundary  in  order  to  illustrate  our  choice  of 
force  density.  Consider  a  boundary  made  up  of  a  continuum  of  material  points  X[s,t)  = 
(x(s,l),y(s,t))  which  is  in  equilibrium  when  its  material  points  coincide  with  the  tether 
points  X*(s)  =  (z'(s),t/*(s)).  We  imagine  there  is  a  spring  connecting  X(s,t)  to  X*(s)  and 
that  it  has  resting  length  zero. 

The  elastic  force  acting  on  a  point  X(s,  t)  is  then: 
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f  (s,t)ds  =  - S[X(s,t )  -  X*(s)]ds 


(4) 


S  is  the  stiffness  constant  of  the  spring,  and  it  may  be  thought  of  in  two  ways: 

1)  A  physical  parameter  which  reflects  material  properties  of  the  elastic  boundary. 

2)  A  numerical  parameter  which  should  be  made  as  large  as  possible  to  strictly  enforce 
the  equilibrium  position  of  the  elastic  boundary. 

This  representation  of  the  force  density  assumes  that  we  know  a  priori  the  equilibrium 
position  of  the  boundary  in  fixed  space.  This  will  be  useful  in  our  test  examples  where  our 
immersed  boundary  will  be  either  a  fixed  cylinder  or  a  fixed  flat  plate.  However,  in  most 
applications,  the  motion  which  we  choose  to  specify  is  the  time-dependent  configuration  of 
the  boundary  with  respect  to  itself.  The  actual  displacement  relative  to  spatial  coordinates 
is  not  preset,  but  is  determined  by  the  equations  of  motion.  For  instance,  we  can  choose 
forces  at  a  boundary  point  due  to  the  rest  of  the  boundary  as: 

1)  Spring-like  forces  which  asks  that  the  "links”  between  successive  points  resist  com¬ 
pression  or  expansion  from  a  given  arclength. 

2)  Bending-resistant  forces  which  asks  that  the  angle  formed  by  neighboring  links  be  a 
given  function  which  changes  with  position  and  time. 

We  refer  the  reader  to  [9]  [16]  for  more  details  on  how  different  boundary  configurations 
can  be  enforced  by  the  proper  choice  of  force  density  functions. 
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3.  Vortex  Method 

In  this  section,  we  will  outline  Chorin’s  vortex  method  [6],  a  grid-free  scheme  which 
represents  the  fluid  field  by  a  discrete  collection  of  vortex  elements. 

Inviscid,  incompressible  flow  in  two  dimensions  in  the  absence  of  external  forces  is  gov¬ 
erned  by  the  Euler  equations: 


Du  „ 

Dt  = 

(5) 

V-u  =  0 

(6) 

The  first  term  in  equation  (5)  is  the  total  derivative.  We  introduce  the  vorticity  i o  = 
V  x  u  •  (0, 0, 1)  =  vx  —  uy  .  Taking  curl  of  equation  (5)  gives: 


(7) 


The  vorticity  acts  as  "particles”  which  move  with  the  fluid.  Since  this  flow  is  inviscid, 
there  is  no  mechanism  for  creation  or  destruction  of  vorticity  (Kelvin’s  circulation  theorem). 
The  vorticity  is  represented  by  a  sum  of  m  point  vortices,  the  i  th  one  centered  at  point  x,- 


with  strength 


m  m 

^  =  J2uf{  =  ^2KiS(x~  x«)  (8) 

i=l  i=  1 

Here  8  is  the  Dirac  delta  function. 

The  incompressiL'lity  condition  implies  the  existence  of  a  stream  function  ip(x,  y)  such 
that: 


(9) 
(10) 
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%  =  u 
ipx  =  -v 


The  relation  between  the  stream  function  and  the  vorticity  is  therefore: 


A  tJ)  =  -w  (11) 

Using  the  representation  (8)  of  vorticity,  equation  (11)  can  be  solved  explicitly  (in  the  whole 
plane)  by  the  formula: 


4>(x)  =  £  ^log\x  ~  x‘ 


The  velocity  can  be  calculated  directly  from  the  stream  function: 


(12) 


«*(x)  =  £ 


(~(P  ~  y.),  (x  ~  ^)) 


«=i 


2tt 


|x-x,| 


(13) 


This  formula,  unfortunately,  is  infinite  at  each  of  the  vortex  centers  x,-  .  This  is  over¬ 


come  by  the  convention  that  the  induced  self-velocity  of  a  point  vortex,  in  the  absence  of 
boundaries,  is  zero.  The  convection  equation  for  vorticity  (7)  is  solved  by  moving  the  point 
vortices  at  the  local  fluid  velocity: 


Ki  (  (y»  Ui)t  ixi  ~  xj)) 

dt  i+i'2-*  lx«'-xil2  1 

Notice  that  this  Lagrangian  description  of  the  vorticity  field  overcomes  the  difficulty  of 

the  nonlinearity  in  the  convection  equation. 

Thus  far  we’ve  been  considering  an  unbounded  domain.  In  a  region  of  fluid  with  sta¬ 
tionary  boundary  dQ,  the  boundary  condition  which  must  be  satisfied  in  this  inviscid  flow 
is: 


u  •  n  =  0  on  dQ.  (15) 

wheie  n  is  the  normal  to  the  boundary.  (If  the  boundary  is  moving  at  velocity  V,  the 
condition  becomes  u  •  n  =  V  ■  n  .)  This  can  be  satisfied  by  adding  a  potential  flow  to  the 
velocity  in  (13): 
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u  =  U  u,  +  V<$ 


(16) 


Taking  divergence  of  this  expression,  we  get: 


A<?  =  0  (17) 

dn<f>  =  — uw  •  n  (18) 

A  Neumann  problem  must  be  solved  for  (f>.  Depending  upon  the  geometry  of  the  problem, 
this  can  be  done  through  the  use  of  conformal  mapping  where  proper  image  vortices  are 
introduced.  [11][19] 

The  vortex  method  for  incompressible,  viscous  flow  is  based  on  these  ideas.  The  major 
difference  between  the  Navier-Stokes  equations  and  the  Euler  equations,  besides  the  diffusion 
term,  is  the  role  of  boundaries.  The  condition  at  a  stationary  boundary  is  now  the  no-slip 
condition: 


u  =  0  on  dQ.  (19) 

Not  only  is  the  normal  velocity  zero,  but  also  the  tangential  velocity.  Also,  unlike  inviscid 
flow,  vorticity  can  be  created. 

Consider  flow  past  an  infinite  flat  plate  with  uniform  velocity  (17,0)  at  infinity.  The 
solution  of  the  inviscid  equations  is  u  =  (£/,  0)  everywhere.  However,  in  the  viscous  case, 
there  must  be  zero  velocity  at  the  plate. 

The  transition  region  from  the  boundary  to  the  uniform  flow  at  infinity  (the  boundary 
layer)  has  thickness  on  the  order  of  1  />//?,  where  R  =  Reynolds  number.  Within  this 
boundary  layer,  vorticity  is  created.  For  large  Reynolds  numbers,  the  boundary  layer  is 
quite  small.  Here,  again,  the  inadequacy  of  finite  difference  schemes  for  large  Reynolds 
number  flow  becomes  evident.  In  order  to  resolve  what  is  going  on  in  the  boundary  layer,  a 
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minimum  number  of  mesh  points  is  needed  within  that  region.  This  leads  to  a  prohibitively 
small  mesh  width  and  time  step. 

Vortex  Method  for  Viscous  flow 

The  Navier-Stokes  equations  for  two-dimensional  flow  in  stream  function  -  vorticity  for¬ 
mulation  are  : 


Dlj 

1 

(20) 

~Dt 

=  nAu 

R 

=  — w 

(21) 

u 

= 

(22) 

V 

=  — 

(23) 

As  before,  the  vorticity  field  is  represented  by  a  discrete  collection  of  vortex  elements. 
In  practice,  it  is  necessary  to  spread  the  point  vortices  into  *  vortex  blobs”  [11].  We  use  the 
approximation  suggested  by  Krasny  [13].  The  velocity  field  in  an  unbounded  fluid  induced 
by  these  modified  point  vortices  would  be: 

|x-X,f  +  <T’  1  ’ 

where  a  is  smoothing  parameter. 

The  no-slip  boundary  condition  on  89.  is  split  up  into  two  parts: 


u  -  n  =  0  (25) 

u-r  =  0  (26) 

Here  n  is  the  normal  vector  and  r  is  the  tangent  vector  to  89. .  The  first  condition  is  satisfied 
by  the  introduction  of  a  potential  flow,  as  explained  above.  The  second  condition  is  satisfied 
by  the  creation  of  vortex  elements  at  the  boundary  with  strength  — u  -  r  per  unit  length  of 
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the  boundary.  (The  boundary  is  partitioned  into  pieces  of  length  As,  and  a  modified  point 
vovtex  of  strength  — 2(u-r)As  is  placed  at  its  center.  The  factor  of  two  accounts  for  elements 
which  will  immediately  fall  outside  the  boundary  during  a  random  walk,  as  explained  below.) 
This  numerical  creation  of  vorticity  at  the  boundary  mimics  what  physically  occurs. 

The  vortex  algorithm  for  solution  of  the  Navier-Stokes  equations  is  as  follows.  Given  the 
centers  x,  and  strengths  «,•  of  the  vortex  elements  at  time  step  n  : 

1.  Satisfy  the  tangential  boundary  condition  by  creating  new  vortex  elements  at  the 
boundary. 

2.  Move  all  of  the  vortex  elements: 

x?+1  =x?  +  A  t(vC  +  V<fin)  +  m  .  (27) 

Here  the  tj,  are  independent  Gaussian  random  variables  with  mean  zero  and  variance  2A t/R. 
This  random  walk  accounts  for  the  diffusion  term  in  the  vorticity  convection  equation. 
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4.  Grid-free  method  with  immersed  boundary 

The  vortex  method  described  above  is  grid-free,  and  thus  the  small  mesh  width  restriction 
of  finite  difference  schemes  for  large  Reynolds  numbers  is  circumvented.  We  have  developed 
a  vortex  method  to  solve  the  coupled  system  of  a  viscous  fluid  and  an  immersed  elastic 
boundary.  Our  treatment  of  the  elastic  boundary  as  a  vorticity  source  is  based  on  McCracken 
and  Peskin  [15]. 

The  fluid  only  feels  the  presence  of  the  immersed  boundary  as  a  singular  distribution  of 
external  forces: 


F (x,<)  =  /  f(s,t)6(x  -  X(s,t))ds  ,  (28) 

J  B 

Here  the  integration  is  over  the  immersed  boundary,  and  6  is  the  two-dimensional  delta 
function.  We  need  to  examine  how  this  contributes  to  the  evolution  of  the  vorticity.  Following 
McCracken  and  Peskin  [15],  we  see  that  the  curl  of  the  external  force  acts  as  a  source  of 
vorticity: 


Doj  1  .  , 

Dt  =RAu  +  (VxF>'k 

(29) 

Here  k  = 

(0,0,1).  We  have: 

(V  x  F)  •  k  =  J  [V  x  (S(x  -  X(s,t))  f)  •  k  is 

(30) 

(31) 

where  f  = 

(/i  j  fi)1  The  above  integrand  can  be  interpreted  as  a  directional  derivative.  We 

discretize  this  integral  as  follows: 

v-«(*+  -  X(s,t))  -  S(x-  -  X(s, 

£  2  h  |f|As 

(32) 

where 

x+  =  X(s,t)  + 

(33) 

x-  =  X(M)  - 

(34) 
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Here  h  Is  a  small  parameter. 

Therefore,  the  immersed  boundary  acts  a  source  of  vortex  dipoles  with  dipole  moment: 

(/2,-/i)AsA<  (35) 

We  multiply  by  At  to  get  the  amount  of  vorticity  to  be  generated  over  one  time  step. 

This  dipole  moment  is  not,  in  general,  normal  to  the  immersed  boundary.  McCracken 
and  Peskin  [15]  split  up  the  force  density  f  into  its  tangential  and  normal  components: 

f  =  frT  +  ft)T)  (36) 

where  r  is  the  unit  tangent  to  the  immersed  boundary  curve,  and  r?  is  the  unit  normal.  It  can 
be  shown,  using  integration  by  parts,  that  the  normal  component  contributes  a  monopole 
layer  of  vorticity  along  the  boundary  with  strength  density: 

^  r  f*)(s)  I 

dsl|X'(s)|J 

The  tangential  component  contributes  a  dipole  layer  as  described  above.  The  process  of 
splitting  up  the  creation  of  vortex  elements  into  a  monopole  layer  and  dipole  layer  oriented 
normal  to  the  boundary  adds  more  numerical  stability  to  the  calculations. 

We  now  outline  our  basic  vortex  method  in  the  plane  within  which  there  is  an  immersed 
boundary.  The  state  of  the  system  at  time  t  =  nAt  is  given  by  the  configuration  X£ 
of  material  points  of  the  immersed  boundary,  and  the  discrete  collection  of  modified  point 
vortices  centered  at  x"  with  strengths  n1-.  The  basic  algorithm  is: 

1.  Compute  force  density  f£  using  the  current  boundary  configuration. 

2.  Use  (1)  to  create  new  vortex  elements  at  the  boundary. 

3.  Move  each  of  the  vortex  elements  in  the  velocity  field  induced  by  the  other  vortex 
elements  plus  a  random  walk  to  simulate  diffusion. 

4.  Compute  the  velocity  at  each  of  the  points  of  the  immersed  boundary  using  (24). 
Move  each  of  these  at  this  local  fluid  velocity. 
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This  algorithm  has  the  advantage  of  being  grid-free,  and  since  it  has  no  real  boundaries, 
we  needn’t  worry  about  the  introduction  of  image  vortices  to  satisfy  the  normal  boundary 
condition  discussed  in  the  previous  section. 

Explicit  evaluation  of' the  force  density  at  the  current  boundary  configuration  produces 
large  instabilities.  The  reason  for  this  is  that  the  boundary  is  stiff  and  the  forces  are 
large  enough  to  make  the  boundary  overshoot  equilibrium  in  one  time  step.  We  use  the 
approximate-implicit  method  described  by  Peskin  [16].  In  general  when  boundary  points 
are  coupled  to  each  other  through  the  forces,  this  approximate-implicit  method  requires  the 
solution  of  a  nonlinear  optimization  problem.  However,  in  the  examples  addressed  in  this 
report,  our  immersed  boundary  points  act  as  independent  springs.  Therefore,  step  (1)  is 
performed  very  quickly. 

In  step  (2),  we  either  create  two  or  three  vortex  elements  per  boundary  point.  This 
depends  on  whether  we  wish  the  dipole  layer  to  be  normal  to  the  boundary  -  in  which  case 
we  create  both  a  dipole  and  a  monopole  layer  of  vorticity. 

Steps  (3)  and  (4)  require  the  computation  of  velocities  of  the  vortex  elements  and  bound¬ 
ary  points  respectively  using  the  current  positions  of  vortex  elements  and  their  strengths. 
We  use  a  second  order  Runge  Kutta  method  to  update  the  positions  of  the  vortex  elements 
and  the  boundary  points. 
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5.  Computational  Results 

Flow  past  a  cylinder 

Our  first  test  problem  will  be  flow  past  a  circular  cylinder.  The  cylinder  is  modelled 
by  a  collection  of  immersed  boundary  points  which  are  tethered  to  fixed  spatial  positions 
by  springs  of  resting  length  zero.  We  are  deliberatly  choosing  a  test  problem  where  the 
equilibrium  position  of  the  immersed  boundary  is  known  and  is  not  time  dependent.  Flow 
past  a  cylinder  has  been  well-studied  since  the  qualitative  features  of  the  flow  (such  as  vortex 
shedding)  are  Reynolds  number  dependent.  For  a  detailed  discussion  of  the  simulation  of 
flow  past  a  cylinder  using  vortex  methods,  the  reader  is  referred  to  Cheer  [4],  Our  simulations 
differ  in  that  the  cylinder  is  represented  as  a  singular  external  force  field  in  the  fluid,  and 
not  a  computational  boundary. 

We  found  that  introducing  the  dipole  layer  normal  to  the  cylinder  made  the  calculations 
more  stable.  Therefore,  at  each  time  step,  we  create  three  vortex  elements  per  boundary 
point.  This  is  shown  in  Figure  1  for  an  immersed  cylinder  made  up  of  forty  immersed 
boundary  points.  The  cylinder  has  a  diameter  d  =  .1  cm  and  it  placed  in  a  uniform  flow 
U  =  2  cm/ sec.  We  used  a  time  step  At  =  .00025  sec,  and  spacing  between  boundary  points 
of  As  =  .Itt/40  cm.  The  numerical  parameter  which  represents  the  distance  of  the  dipole 
vortex  elements  from  the  immersed  boundary  is  h  =  .8As.  The  viscosity  was  chosen  to  be 
10“8.  The  smoothing  parameter  a  =  .75As.  The  Reynolds  number  based  upon  the  uniform 
flow  velocity  and  diameter  of  the  cylinder  is  therefore  on  the  order  of  107.  The  tethering 
stiffness  constant  used  was  106. 

Figures  2  and  3  show  velocity  fields  of  the  flow  about  the  cylinder  at  time  t  =  .0125sec. 
We  show  two  space  scales  to  emphasize  the  fact  that  because  this  is  a  grid-free  method,  the 
velocities  may  be  calculated  from  the  vortex  element  configuration  at  any  spatial  points  by 
using  (24).  We  are  therefore  able  to  resolve  the  boundary  layer  with  a  much  finer  degree  of 
accuracy  than  if  we  were  working  on  a  fixed  grid.  Smoothing  does  take  place,  however,  both 
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in  the  representation  of  the  vortex  blob  functions  (parameter  a)  and  in  the  finiteness  of  the 
dipole  layer  (parameter  h ). 

Stagnation  points  at  the  upstream  and  downstream  centerline  of  the  cylinder  can  be  seen 
in  Figures  2  and  3.  The  magnitude  of  the  velocity  vectors  at  the  top  and  bottom  of  the 
cylinder  are  larger  than  the  free  stream  velocity,  as  is  to  be  expected.  The  flowfield  within 
the  cylinder  is  nearly  zero,  which  is  evidence  that  the  tethering  forces  are  doing  their  job. 
However,  these  flowfields  do  not  yet  show  vortex  formation  and  shedding  at  the  cylinder. 
We  have  not  run  this  code  long  enough  to  exhibit  these  features,  which  become  evident  at 
such  a  high  Reynolds  number.  In  order  for  the  approximate-implicit  force  calculations  to 
remain  stable,  we  need  to  use  a  relatively  small  time  step.  This  causes  the  creation  of  an 
extremely  large  number  of  vortex  elements  in  a  small  real-time  simulation.  We  will  discuss 
some  remedies  to  this  in  Section  6. 

Flow  past  a  flat  plate 

In  this  example,  we  will  simulate  incompressible,  inviscid  flow  about  a  flat  plate  which 
is  oriented  normal  to  the  free  flow.  There  are  two  possible  inviscid  flows  which  are  solu¬ 
tions.  One  is  a  continuous  potential  flow  which  is  steady  and  symmetric  with  respect  to  the 
plate.  However,  the  velocities  near  the  ends  of  the  plate  are  infinite.  Another  solution  is 
an  unsteady,  discontinuous  potential  flow,  which  sheds  vorticity  from  the  plate.  Here,  the 
velocities  at  the  ends  of  the  plate  are  finite.  This  solution,  when  compared  with  experiments, 
is  more  realistic  and  can  be  thought  of  as  the  vanishing  viscosity  limit.  [14] 

The  flat  plate  is  represented  by  an  elastic  boundary  which  is  made  up  of  forty  immersed 
boundary  points.  The  plate  lies  on  the  line  x  =  1  and  goes  from  y  =  1  to  y  =  3.  It  is 
tethered  to  space  with  springs  of  resting  length  zero  and  stiffness  constants  107.  A  uniform 
background  flow  of  U  =  .5  is  imposed.  The  numerical  parameters  used  are  As  =  .05, 
h  —  .4As,  and  a  =  .75As. 

To  modify  the  algorithm  in  order  to  model  invisid  flow,  we  discard  the  random  walk 
which  approximates  diffusion.  Moreover,  we  only  force  the  normal  component  of  velocity  of 
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the  plate  to  coincide  with  that  of  the  fluid.  (Our  immersed  boundary  point  can  therefore 
only  move  in  the  x—  direction.)  In  inviscid  flow,  tangential  slip  is  allowed. 

In  these  simulations,  we  are  only  going  to  create  a  single  dipole  layer  of  vorticity  at  the 
plate  ,i.e.  two  vortex  elements  per  boundary  point  are  created  at  each  time  step.  This  dipole 
layer  is  actually  tangent  to  the  plate. 

Figure  4  shows  velocity  fields  about  the  plate  at  t  =  .04,  <  =  .25,  t  =  .5,  and  t  =  .75.  The 
first  flow  field  appears  to  be  symmetric  about  the  plate;  a  wake  has  not  yet  had  a  chance  to 
form.  The  successive  flow  fields  show  the  formation  of  two  counter-rotating  vortices  at  the 
ends  of  the  plate.  These  are  getting  bigger  as  time  goes  on. 

Our  numerical  scheme  has  picked  out  the  unsteady,  inviscid  solution.  We  believe  this 
happens  because  our  velocity  field,  due  to  the  smoothing  parameter  a,  is  not  allowed  to  be 
infinite  anywhere.  Moreover,  the  finiteness  of  our  dipole  layer  h  adds  some  smearing  which 
is  manifested  as  numerical  viscosity. 

The  dipole  layer  of  vortex  elements  is  initialized  to  lie  tangentially  along  the  plate,  except 
for  the  two  elements  which  are  placed  a  distance  of  h  from  the  top  and  the  bottom.  The 
vortex  elements  on  the  plate  move  approximately  with  the  same  velocities  as  the  tethered 
boundary  points,  and  thus,  approximately  remain  on  the  plate.  However,  vortex  elements 
are  shed  from  the  top  and  bottom.  Figure  5  shows  the  location  of  vortex  elements  at  time 
t  =  .04  and  at  t  =  .25.  Figure  6  shows  a  contour  plot  of  vorticitj  at  t  —  .25.  These  figures 
qualitatively  compare  tc  the  vortex  sheet  roll-up  described  in  Krasny  [14]. 


6.  Conclusions  and  future  work 

Finite  difference  schemes  do  not  perform  well  in  simulations  of  high  Reynolds  number 
flow  because  a  restrictive  u».  .*  pr  of  grid  points  must  be  used  to  resolve  boundary  layers.  We 
have  presented  a  grid-f,  *  *va.  ical  method  which  may  be  used  to  calculate  flows  around  an 
immersed  elastic  bourn.  *y.  ’•> ;  nave  presented  numerical  results  in  two  simple  cases,  where 
the  boundary  motion  is  sp  cified  and  is  not  time  dependent.  Our  initial  results  show  that 
the  method  gives  the  com  ,  qualitative  features  of  the  desired  flow. 

The  major  drawback  o i  '.ms  algorithm  is  the  growth  of  the  number  of  vortex  elements  at 
each  time  step.  The  computation  of  the  vortex- vortex  interactions  becomes  very  expensive. 
Presently,  we  are  performing  the  straightforward  0{m\)  operations,  where  mn  is  the  total 
number  of  vortex  elements  present  at  time  step  n.  There  are,  however,  techniques  which 
can  be  used  to  speed  up  this  calculation:  Anderson’s  method  of  local  corrections  [1],  or  the 
multipole  method  of  Greengard  and  Rokhlin  [12] .  We  are  planing  to  incorporate  one  of  these 
in  our  code. 

As  the  algorithm  stands,  mn  is  always  increasing  by  two  or  three  times  the  number  of 
boundary  p-.uits  each  time  step.  In  an  exterior  problem,  like  flow  past  a  fiat  plate,  we  can 
throw  away  vortex  elements  when  they  move  a  certain  distance  from  the  immersed  boundary. 
However,  in  interior  calculations,  like  blood  flow  in  the  heart,  some  vortex  elements  are 
trapped  inside  the  immersed  boundary.  We  need  to  investigate  the  systematic  merging  of 
nearby  vortex  elements  in  order  to  make  large  time  simulations  feasible. 

We  need  to  establish  the  dependence  of  this  numerical  method  on  the  numerical  param¬ 
eters.  These  are: 

1)  The  equilibrium  distance  between  boundary  points  As. 

2)  The  small  parameter  h  which  determines  the  distance  of  the  dipole  layer  from  the 
•boundary. 

3)  The  parameter  cr  used  in  smoothing  the  velocity  held  near  the  center  of  the  vortex 
blobs. 
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We  would  like  to  show  convergence  of  this  numerical  scheme  as  these  parameters  are 
refined.  We  should  also  include  another  parameter,  Kmax  which  would  be  the  maximum 
strength  any  one  vortex  element  can  carry.  If  a  strength- >  Kmax  is  calculated,  that  vortex 
element  ought  to  be  split  up  into  other  elements  of  smaller  strengths,  which  add  up  to  /c,-. 
This  has  shown  to  be  very  important  in  the  convergence  of  standard  vortex  methods  [21]. 

The  numerical  examples  presented  here  were  chosen  for  their  simplicity.  These  immersed 
boundaries  do  not  undergo  time  dependent  mol  .ons,  and  their  equilibri’im  position  in  space 
has  been  specified  by  specifying  tether  points.  Moreover,  each  boundary  point  is  uncou¬ 
pled  from  its  neighbors  in  the  force  density  calculations.  Once  a  better  understanding  of 
the  numerical  method  is  achieved  on  these  test  problems,  we  can  use  it,  with  very  minor 
modifications,  to  model  time  dependent  problems  in  biofluiddynamics. 
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FLOW  PAST  FLAT  PLATE 


POSITION  OF  VORTEX  ELEMENTS 


POSITION  OF  VORTEX  ELEMENTS 


FIGURE  5 
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